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We consider a protocol to perform the optimal quantum state discrimination of N linearly in- 
dependent non-orthogonal pure quantum states and present a computational code. Through the 
extension of the original Hilbert space, it is possible to perform an unitary operation yielding a final 
configuration, which gives the best discrimination without ambiguity by means of von Neumann 
measurements. Our goal is to introduce a detailed general mathematical procedure to realize this 
task by means of semidefinite programming and norm minimization. The former is used to fix which 
is the best detection probability amplitude for each state of the ensemble. The latter determines the 
matrix which leads the states to the final configuration. In a final step, we decompose the unitary 
transformation in a sequence of two-level rotation matrices. 

PACS numbers: 



I. INTRODUCTION 

Quantum state discrimination is a crucial prob- 
lem in quantum information theory, and it is well 
known the impossibility of doing this perfectly for 
non-orthogonal states [Jj. Nevertheless, if a non zero 
probability of inconclusive results is allowed for, it 
is possible to never mistake a state for another, by 
means of an appropriate Positive Operator Valued 
Measure (POVM). This strategy is known as unam- 
biguous state discrimination (USD), and the best pro- 
cedure of this kind is that which minimizes the prob- 
ability of inconclusive results. The discrimination 
of two equally probable non-orthogonal pure states 
was firstly considered by Ivanovich Q, Dieks Q and 
Peres Q . The case of two states with unequal prior 
probabilities was treated by Jaeger and Shimony 0. 
Chefles showed that USD of N pure quantum 
states is possible if and only if they are linearly in- 
dependent. 

With the aid of Semidefinite Programming (SDP), 
Eldar Q obtained the set of necessary and sufficient 
conditions the USD measurement operators satisfy. 
She also showed that for a given set of pure states, 
there always is an ensemble for which the optimal 
USD detects each state with equal probability (i.e, 
an Equal Probability Measurement - EPM). 

In references [TIH and EL 

the USDproblem is ap- 
proached via the Neumark theorem [y, i.e, the real- 
ization of a POVM by means of projective measure- 
ments in an extended Hilbert space. In it is done 
in the context of linear quantum optics , whereas in 
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(Tll | an ions trap architecture is considered. 

Our main goal is to derive the transformation which 
maps N non-orthogonal pure states in a set of states 
that can be discriminated by usual projective mea- 
surements in an extended Hilbert space. It is equiv- 
alent to a generalized measurement in the original 
space, and is the content of Neumark's theorem. We 
then present a simple MATLAB code that takes as 
input the ensemble of non-orthogonal states and out- 
puts the best set of discriminable states and the perti- 
nent transformations. We believe this computational 
code can be pedagogically useful to anyone getting 
started with quantum information and semidefinite 
programming. 

The paper is organized as follows. In section I 
we discuss generalized measurements and semidefinite 
programming. In section II we present a protocol to 
discriminate pure non-orthogonal states and, in sec- 
tion IV, we exercise with a numerical example in the 
context of quantum key distribution. We conclude in 
section V. In the appendix, we furnish our MATLAB 
code. 



II. POVMS AND DISCRIMINATION WITH 
SDP 

The operation commonly known as positive 
operator-valued measure can be performed with a set 
of quantum detection operators lik, where the prob- 
ability pt of obtaining the state labeled by the index 
k is given by 

p k = Tr(n k p), (1) 

where p is the density operator of the system. As 
the probabilities are obviously non-negative reals and 
sum to one, all the quantum detection operators are 
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semidefinite operators and form a resolution of the 
identity, 

E n ^ = / - ( 2 ) 

k 

The problem of USD for N pure non-orthogonal states 
can be stated as follows. We assume that a quan- 
tum system was prepared in one of the N pure states 
{|Qi)}j i — 1,.-,N. Each state |Q,) is in a N- 
dimensional Hilbert space. In order to identify a 
state or to return an inconclusive result, the measure- 
ment operators must obey (Qj| Ilfc \Qi) = PiSik with 
< Pi < 1. Therefore, to satisfy condition (2J, we 
have Yli=o ^» = ^ an< ^ the inconclusive result is given 
by IIo=/-££ini. 

In [8], Eldar showed that the optimal unambiguous 
discrimination can be formulated as a SDP problem. 
The measurement operators are expressed in the form 

ni=pi\Qi){Qi\, (3) 

where each state |Qj) is in a iV-dimensional Hilbert 
space, and the IL are not orthogonal projectors in 
this space. The vectors |Qj) are the reciprocal states 
associated with \Qi), such that 

(QilQfc) = QkSik , l<i,k<N (4) 

where the factor Qk indicates the scalar product is 
not normalized. Given the matrix 'J, whose columns 
are the vectors \Qi), the \Qi) are the columns of the 
matrix 'J, namely, 

§ = (5) 

Given an ensemble p, where each state \Qi) is pre- 
pared with probability fn, the total probability of a 
successful detection is 

N N 

p d = ^/i i (Q,|n l |Q l ) = ^inpi. (6) 

i=l i=l 

The problem of optimal USD is then to find measure- 
ment operators H, or equivalently, the probabilities 
Pi , which maximize Pd, subject to the constraint {pj, 
which can be recast as 

N 

I-^JH\Qi){Qi\>0- (7) 

i=l 

A SDP problem is to find x € R m , which minimizes 
the linear function c T x , subject to the matrix in- 
equality F(x) = Fq + J2iLi XiFi > 0, where the prob- 
lem data are the vector c £ R m and the m+1 complex 



Hermitian matrices Fi jl2j. It is known as the primal 
formulation of SDP. 

Equations © and {JJ can be recast as a SDP prob- 
lem, namely: 

min peJi N{-/i T p}, (8) 
subject to N + 1 constraints, 

AT 

/-X>I&X&I >o, 

i=l 

Pi > 0, l<i<N. (9) 

III. PROTOCOL FOR OPTIMAL 
DISCRIMINATION OF PURE STATES VIA 
ROTATIONS AND VON NEUMANN 
MEASUREMENTS 

We start out by rewriting the N entry states to 
be discriminated in a ladder form in the orthonormal 
basis {\i) , i = 1, . . . , N}, 

\Q±) = |i>, 

|Q 2 ) = c 12 |l)+c 22 |2), 

\Qs) = Ci 3 |l)+C 23 |2)+ C 33|3) 



\Qn) = ci^|1) + ... + cjw|JV) , (10) 

where € Hn (TV-dimensional Hilbert space). 

This can always be done by means of a unitary trans- 
formation Uq. Then, to apply Neumark's theorem, we 
extend the original Hilbert space, e.g. through addi- 
tion of ancillas, to 2N — 1 dimensions and map the 
original states to the final configuration, 

\Qif) = 9x\l)+9mi\N+l)+...+g 2 ^-i\2N-l) 
\Qif) = 92 |2) +g 2N \N+l)+... +g 3N _ 2 \2N-l) 
\Q»f) = 93 1 3) +g 3 N-i\N+l)+...+g m ^ 2 N- 2) 

\Qif) = 9i\i) + 9[±(2N+3^)^i]\ N + !> + 

• • • + ff[JV(i+l)+i(l-i)-l] \ 2N + 1 - i) 

\QNj) = .9Ar|A0+.g [ i A r ( ^ 3) _ 1] |A r + l>. (11) 

Now, a projective measurement in the orthonormal 
basis {\i) ,i = 1, . . . , 2N — 1} yields an unambiguous 
discrimination. The state labeled by i is identified 
with probability gf , when the measurement collapses 
to \i) for 1 < i < N. For other values of i the result is 
inconclusive. Therefore, the first TV gi are chosen to 
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produce the best possible discrimination, i.e, gi = y/pl 
for the pi defined in Eq.3, which are determined by 
SDP. The other coefficients {g it i = N + 1 to ~(N + 
3) — 1} are fixed in order to preserve normalization 
and the scalar products among the original states. 

Once the final configuration is known, we want to 
determine the unitary transformation U± which maps 
the original states to it. Norm minimization 0| im- 
plies, 



N 

E 



\Qif) -Ui\Qi)\\ 



where 



2N - 2Re 



N 



Tr(A'ut)] = 0, (12) 



(13) 



A' is a singular (2N—1) x (2N—1) complex matrix. 
The singular value decomposition (SVD) of A' is 
A' = VYW\ with V and W (2JV-1) x (2N-1) unitary 
matrices, and S a diagonal {2N— l)x (2N— 1) matrix, 
then Ui is given by 



Ui = vwl 



(14) 



The unitary transformation U = UxUq can now be 
decomposed [15] in a sequence of rotations (Rkl) in 
the hyperplanes (kl), as 



N-i N 



U 



n n R i 

k=l l=k+l 



(15) 



The Rki are two-level matrices, with the four non triv- 
ial entries, 



[R 



kllkk 



PYkk 



[Rkl]k 



[Rki]u - -\Rki\kk-, 
[U]ik 



[Rki]ik - [Rki]*ki- 

The trivial entries are 

[Rkl]mm = 1,™ / [k,l), 



[Rk 



0,(m,n)^(k,l). 



We can also express the rotations R^ kl in terms of 
Pauli matrices, i£j(#) = exp(— i0a z /2) (the super- 
index k indicates the rotation is over the ket |fc)), 
R l y (9) = exp(-i9a y /2). Hence, 

Rl(a,/3, 7,6) = e ia R k MR l y{l)R k M. (16) 



Therefore the unitary transformation U = U\Uq, 
followed by a projective measurement in the basis 
{\i),i = 1, ...,2N — 1}, discriminates unambiguously 
any N pure non-orthogonal states. Summarizing, we 
have the following algorithm: 

• Rewrite the entry states in a ladder form (U ). 

• Fix the conclusive amplitudes (gi to gjy) using 
SDP. 

• Fix the inconclusive amplitudes (gN+i to 
<?«(at+3)_i) such that normalization and scalar 
products among the states be preserved. 

• Build the unitary transformation (U\) that 
maps the entry states in the ladder form to the 
final discriminable configuration. 

• Decompose U = U\ Uq as a product of one-qudit 
rotations. 



IV. EXAMPLE 

As an example of application of our code, we con- 
sider how a spy (Eve) could use USD to eavesdrop 
E two parties (Alice and Bob) establishing a 
BB84 cryptographic key. 

Alice and Bob are communicating by means of low 
intensity laser pulses in a network of lossy optical 
fibers. The signals sent by Alice are coherent states 
with a certain mean number of photons (fx). The 
polarization of the pulses are randomly chosen with 
equal probability among four possibilities, namely, di- 
agonal to right or left (\d+) , \d— )), and circularly 
polarized to right or left (|c— ) , |c+)). Bob directs 
the pulses he receives to detectors, which do not dis- 
tinguish photon numbers, preceded by polarization 
analyzers which can be either for linear or circular 
polarization, also chosen randomly with equal prob- 
ability. After a certain number of pulses, Alice pub- 
licly announces the sequence of polarization basis she 
used. Bob then checks which pulses he detected using 
the compatible polarization analyzer. After discount- 
ing dark counts and losses in the fibers, the matching 
detections allow Alice and Bob to establish a crypto- 
graphic key, i.e., a long sequence of zeroes and ones 
(say \d+), \c+) for one and \d—), |c— ) for zero). This 
is the very well known BB84 protocol. 

Eve could try to obtain this key as follows. She 
probes the network and measures the mean number 
of photons in the pulses, which can be done with- 
out disturbing the polarization. Now she knows that 
most of the time Bob receives states of the type 
\polarization)®^ . Suppose fi is 3. Eve then pre- 
pares an USD scheme for the states |d+)|d+)|d+)j 
\d-)\d-)\d-), \c+)\c+)\c+), |c-)|c-)|c-). 
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Running our code, all the pertinent parameters to 
prepare the USD are yielded. In particular, we learn 
that Eve correctly identifies the polarization state 
with a probability of 50%, if the state has three pho- 
tons. When she succeeds, she prepares a state with 



the correct polarization and sends it to Bob, who will 
never know it came from Eve. If she fails, she does 
nothing at all, and Bob could think it is a dark count 
or a loss in the network. 
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Let's consider the four quantum states 



\d+)\d+)\d+) 
\d-)\d-)\d-) 



where 



iQi) 

\Q\) 

|Q|) = |c+)|c+)|c+) 
IQl) = | C ->|c->|c-) 



\d+) 
\d-) 



(10) + |1»/V2 
(|0>-|l»A/2 
| C +) = (|0)+z|l))A/2 
|c-> = (|0>-i|l»/V2 . 

Rewrite the entry states in a ladder form 

IQi)=|i) 

IQa> = |2) 

|Q 3 ) = -(0.25- 0.25i)|l) + 

-(0.25 + 0.25i)|2) +0.8660|3) 

|Q 4 ) = -(0.25 + 0.25i)|l) + 

-(0.25- 0.25i)|2) +0.8660|4), 



(17) 
(18) 
(19) 
(20) 



(21) 



where the unitary matrix Uq for the transformation 
is given above. 



The best conclusive probability amplitude given by 
the SDP technique is given by gi=<?2=<73=<74=0.7071, 
where the weight of each state in the ensemble is 
A i i=A i 2=A i 3=A t 4— 0.25. The inconclusive amplitudes 
(35 to 1713) are fixed such that normalization and 
scalar products among the states be preserved. For 
this case, g§ — —0.3536 + 0.3536i, ge — —0.3536 — 
0.3536?, (77 = 0, gs = —0.3536 — 0.3536i, gg = 
-0.3536 + 0.3536i, g w = 0, g n = 0, g 12 = 0.7071 
and gig = 0.7071. 

The final discriminable configuration in the ex- 
tended Hilbert space is 



\Q lf ) = 0.7071|1) - (0.3536 - 0.3536z)|5) 

-(0.3536 + 0.3536i)|6) 
\Q 2f ) = 0.7071|2) - (0.3536 + 0.3536i)|5) 

-(0.3536- 0.3536i)|6) 
\Q 3f ) = 0.7071|3) + 0.7071|6) 
\Qi f )= 0.7071|4) + 0.7071|5) . (22) 



Now we build the unitary transformation (U±) that 
maps the entry states in the ladder form to the final 
discriminable configuration. Following the procedure 
described by equations 11 -ill and (1141 . we obtain: 
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As the last step, we wish to decompose the resultant 
matrix U — U\Uq as a product of one-qudit rotations. 
Therefore, we have 



u=r\ 2 r\ 3 



R 7,8 



where R\ 5 — i?} 6 — i?| 7 — i?} 8 — 4,7— -4,8- 



R. 
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Table I: A specification of the parameters a, /3, 7, and 5 
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The anj 


;les a, /3, 7, and 5 


are given 


in degrees. 




a 


P 


7/2 


8 


Ri 2 


90.0 


0.0 


45.0 


180.0 


Ri,s 


90.0 


0.0 


35.27 


180.0 


Rl,4 


90.0 


0.0 


30.0 


180.0 


#2,3 


-90.0 


-2.2066 


-36.9938 


-128.4034 


R"2 4 


-90.0 


-138.9380 


46.6463 


-41.0619 


Ri 5 


90.0 


160.2880 
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19.7119 


_/?2 6 


-90.0 


-115.9253 


16.3825 


-64.0746 


R 3A 


-90.0 
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R3.5 


-90.0 
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0.0 


32.4564 


180.0 


#4,8 


90.0 
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-"-7,8" 



Specifying the parameters (a, (3, 7, <5) in each step of 
the decomposition we conclude the the unambiguous 
discrimination protocol. The parameter's values are 
given in Table I. 



V. CONCLUSIONS 

We showed a general algorithm to perform the op- 
timal discrimination of TV linearly independent non- 
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orthogonal pure quantum states by means of semidef- 
inite programming and norm minimization. In addi- 
tion, we presented a simple computational code that 
takes as input the ensemble of non-orthogonal states 
and outputs the best set of discriminable states and 
the sequence of two-level rotation matrices. As a nu- 
merical example, we studied an USD attack to the 
BB84 protocol. 
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Appendix A: DETERMINATION OF THE 
INCONCLUSIVE PROBABILITY 
AMPLITUDES 

In this appendix we supply the procedure to cal- 
culate the inconclusive probability amplitudes in the 
final configuration, i.e, gi, for i = N+l to ^N(N+3)— 
1. We detail the procedure for the case N = 3, and 
then discuss the generalization, presenting a MAT- 
LAB code for arbitrary N. 

Let's consider three non-orthogonal quantum states 
in a ladder form. It is convenient to factorize a global 
phase and set the coefficient multiplying the the ket 
|1) as real, for all states, then: 



\Qi) = |i> 
IQa) = oi |i; 
IQs) = a 3 |i; 



a 2 e 



|2> 



a 4 e ie2 |2) +a 5 e l6a |3) , (Al) 



The final discriminable configuration in the extended 
Hilbert space is 



\Qif) = 5i|l>+54|4)+ 55 |5) 
\Q2f) = ff 2|2)+ ff6 |4)+ 57 |5) 
IQs/) = g 3 \3)+g 8 \4). 



(A2) 



The conclusive probability amplitudes 51,32,53 axe 
determined via SDP. The inconclusive probability am- 
plitudes 54 to 58 must be determined preserving nor- 
malization and scalar products among the vectors. 
From normalization we have, 



\9i\ +\9i\ 

M 2 + |.96| 2 + |57| z = I 

|2 1 |„ |2 



\to\ 2 
|2 



133 I 



1^ = 1 



(A3) 
(A4) 
(A5) 



and the scalar products are 

(Q1/IQ2/) = (Qi\Q 2 ), 
(Q1/IQ3/) = (O1IQ3), 

(Q2/IQ3/) = (g 2 |g 3 >. (A6) 

ijAfill can be written in the form 

9496 +9*597 = ai (A7) 
gigs, = a 3 (A8) 
5e38 = a 2 a 4 e l< - 9 ^ + ai a 3 . (A9) 

From l!A8|l it is observed that 54 and 58 can be taken 
as real. Due to ijAfiJ) we have 58 = +yl — \g-i\ 2 - 
Substituting g% in l|A8J) we get 54 = +a 3 /y/l — I53I 2 . 
Substituting 58 in (j A0|) we have 

[a 2 a 4 cos + 0,10,3] — «[a 2 a4 sin 9] 
36 = ■ • (A10) 

v 1 - l33p 

with 6 = 02-01. 

55 and 57 are still undetermined. Writing them as 

35 = gf + igl m e 57 = gf + «37 m , and usin s 1331, 

(EH, (53), then 



5? = ±Jb-(g^y, (All) 



3? 



R R i I m I m 
97 95 + 95 97 = c, 
R Ira Irn R j 

35 37 - 3 5 37 = o, 
with a, 6, c and d reals given by 

b = 1-|32| 2 -|3 6 | 2 



±V a -C95 /m ) 2 , (A12) 
(A13) 
(A14) 



a i 



a3[aia 3 + 0204 cos 0] 



i-N 



d = 



0^0204 sin 
1-|33| 2 



(A15) 



Substituting l|AlT)l and fA12l into (SUt and llA14|l 

we have 

- (5 5 /m ) 2 Jb- (5 7 /m ) 2 + 3 5 /m 37 /m = c , (A16) 

\/«-(35 /m ) 2 3 7 /m - 3s 7m \/&-(37 /m ) 2 = d , (A17) 
From ijAlfiJl we have 

Im _ 2cgi m ± V(2c5i m ) 2 -4(6(5^™)2+ c 2_a6)a _ 



37 



2a 



and from !|A17|) 

{g i m) A + [(4d 2 -2a6)(5|™) 2 + 2ad 2 



+ 



b (gl m ) ~ 6d 2 6(5^ m ) 2 +rf 



-(A18) 

(5 7 /m ) 2 + 
.(A19) 
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Substituting ijlTs)) into tKTfy we find that 
must obey the following polynomial 

A{gi m f + B{gi m f + C7(g 5 /m ) 4 + D(g I 5 m ) 2 + E=0 ,(A20) 

with real coefficients 

A = [{Ad 2 - Aab)(c 2 - ab) - 4c 2 d 2 ] 2 , 

B = [(Ad 2 - Aab){c 2 - ab) + Ac 2 d 2 ] x 

x [8a 2 b(c 2 -ab- d 2 )] - QAc 2 d 2 (c 2 - ab)(2a 2 b) , 

C = [a 2 d 4 + a 2 (c 2 - ab) 2 - 2d 2 a 2 (c 2 - ab)] x 

x [(Ad 2 - Aab)(c 2 - ab) + Ac 2 d 2 ] + 

+ [Aa 2 b(c 2 -ba-d 2 )] 2 + 

+6Ac 2 d 2 (c 2 - ab)(d 2 + ab)a 2 , 

D = [a 2 d 4 + a 2 {c 2 - ab) 2 - 2d 2 a 2 (c 2 - ab)] x 

x [Aa 2 b(c 2 - ba - d 2 )] , 

E = [a 2 d* + a 2 (c 2 - ab) 2 - 2d 2 a 2 (c 2 - a6)] 2 (A21) 

The roots of this polynomial are easily obtained ^3] • 
We choose for any real root, such that < 
(fff™) 2 < a - Once <?5 is determined, we calculate 
g\ m from (JXTHJ. Then llAlT)l and l|AT2")) are used to 
calculate gf and g^. Now all the parameters in the 
final configuration are determined. 

It is straightforward to extend this procedure for 
arbitrary N . Summarizing: we determine the con- 
clusive amplitudes, g\ to gjv, by SDP; then, starting 
from the last inconclusive amplitude, (7[Ajv(jv+3)-i]i 
we use the relations for normalization and scalar prod- 
ucts and thus determine all the remaining parameters, 
except for gw-\ and g^N-2, which are obtained by 
means of the polynomial 

+ D(g{^_ 1 ) 2 +E = Q. (A22) 

7, Unambiguos State Discrimination 

7, Wilson R.M. Rabelo (wilson@unifap.br) 

7. 09/2005 UFMG Quantum Information Group 

7. BEGIN INPUT - BB84 

7oNumber of states 
N=4; 

7,Hilbert space dimension 
7. (at least 2N-1) 
dim=8 ; 

7oProbability of each state in the 
7, ensemble 
Mi=-1*[0.25; 0.25; 0.25; 0.25]; 
zero= [1 ;0] ; 
one= [0 ; 1] ; 

lr= (zero+one) / sqrt (2) ; 



ll=(zero-one)/sqrt(2) ; 
cr= (zero+i*one) / sqrt (2) ; 
cl=(zero-i*one)/sqrt(2) ; 
lr2=kron(lr,lr) ; 

112=kron(ll,ll) 5 
cr2=kron(cr , cr) ; 
cl2=kron(cl , cl) ; 
lr3=kron(lr2,lr) ; 

113=kron(112,ll) 5 
cr3=kron(cr2 , cr) ; 
cl3=kron(cl2 , cl) ; 

7.1nput sates QII= Q_l Q_N 

QII=[lr3 113 cr3 cl3] ; 

7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7. 

7, END INPUT 

7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7. 

7. BEGINNING: 

v/:/:/:/:/:/:/:/:/^^ 

7.N initial states: 
disp('N initial states:') 
disp(QII) 

7.Scalar product for N initial states: 

for jL=l:N 
for iL=l:N 
Esc_QII=QII(l:dim,iL) '*QII(l:dim, jL) ; 
esc_QII(iL, jL)=Esc_QII; 
end 
end 

7. Rewrite the entry states 
7. in ladder form (Uo) 

7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7. 
7.for ig=l:2*N-l 
for ig=l:dim 
for jg=l:N 

c(ig,jg)=0.0; 

end 
end 

c(l,l)=l-0; 
for ittt=l:N-l 
for j=ittt+l:N 
Sf at_w=0.0; 
if ittt > 1 

for xxw=l:ittt-l 
Sf at_w=Sf at_w+conj (c(ittt-xxw, ittt) )*c(ittt-xxw, j ) ; 
end 
end 

if N==2 

c(ittt, j)=esc_QII(ittt, j) ; 
c (ittt+1 , ittt+1) =sqrt (1-conj (c (ittt , ittt+1) ) *c (ittt , itt 
continue 
end 

if N >= 3 

if ittt==l 

c(ittt, j)=esc_QII(ittt, j) ; 
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else nnnnnnnnnnnununnnntn 

c(ittt , j) = (esc_QII(ittt, j)-Sf at_w)/c(ittt ,iffc*rt^jjj = l:N 

end g(jj)=sqrt(P(jj)); 

end end 

end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

Swfat=0.0; disp( 'Conclusive amplitudes via SDP :') 

for xo=l:ittt g 

Swf at=Swf at+conj (c (xo , ittt+1) ) *c (xo , ittt+1) %%%%%%%%%%%%%%%%11%%111%%11%%%11%%%1'm 

end pause 
c(ittt+l,ittt+l)=sqrt(l-Swfat) ; 

end for ii=N+l: ( (N/2) * (N+3) ) -1 

g(ii)=0.0; 

Q=c; end 
format short 

AAA=zeros(dim) ; g(((N/2)*(N+3))-l)=sqrt(l-g(N)~2) ; 
for ux=l:N 

AAA=AAA +Q(l:dim,ux)*QII(l:dim,ux) ' ; x=l; Sf atg=0; j j j=0; 

end w=3 ; 

[V_0,Sigma_0,W_0]=svd(AAA) ; for jL=N:-l:3 

disp('The unitary matrix Uo ') j=((N/2)*(N+3))-x; 

dispCto put states in ladder form :') jjj=jjj + l; 

Uo=V_0*W_0' k=jjj; 

disp('Test Uo, i.e, (Uo*)(Uo)=') f or kL=jL-l : -1 : 1 

disp(Uo*Uo') if kL > 1 

dispC Initial configuration: [Ql. . .QN] ') k=k+l; 

dispC in ladder form: ') j=j -k; 

for iu=l:N if w > 3 

Qcc(: ,iu)=Uo*QII(l:dim,iu) ; Sfatg=0; 

end P=0; 

Q fatg=0; 

pause for f=((w/2)*(w-5))+5: ((w/2)*(w-3))+l 

format long P=P+1; 

raffirmny.rammray.nramy,ram E(f )= g ( ( (n/2) * (n+3) ) -f )/ g ( ( (n/2)* . . . 

"/Scalar product for N states (N+3))-x); 

7, in ladder form: fatg=conj (g(j-P))*E(f ) ; 

for jL=2:N Sf atg=Sf atg+f atg; 

for iL=l:jL-l end 

Esc_Q=Q(l:N,iL) '*Q(l:N,jL) ; end 

esc_Q(iL,jL)=Esc_Q; g(j)=((esc_Q(kL, jL) ) / (g( ( (N/2) * . . . 

end (N+3))-x)))-Sfatg; 

end g(j)=conj (g(j)) ; 

y.r/x/:/x/x/.r/.r/.sDP Ap P roachr/x/x/x/:/.n°/.n ei se 

c=sdpvar(N,l) ; j=j-k; 

p=sdpvar (N, 1) ; if w > 3 

F_0=eye(dim) ; Sfatg=0; 

Qt=Q*inv(Q'*Q) ; P=0; 

D=F_0; fatg=0; 

for ib=l:N for f=((w/2)*(w-5))+5: ((w/2)*(w-3))+l 

D=D-p(ib)*Qt(l:dim,ib)*Qt(l:dim,ib) ' ; P=P+1; 

end ' E(f)=g(((N/2)*(N+3))-f)/g(((N/2)*. . . 

yalmip('info') ; (N+3))-x); 

F= set(D > 0); f atg=conj (g( j -P) ) *E(f ) ; 

for ikp=l:N Sfatg=Sf atg+f atg; 

F=F+set(p(ikp)>0) ; end 

end end 

solvesdp(F,Mi'*p) ; g(j)=((esc_Q(kL, jL))/(g(((N/2)*. . . 

P=double(p) (N+3))-x)))-Sfatg; 
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g(j)=conj (g(j)) ; 
end 
end 
if N==3 

continue 

end 

if jL > 3 
w=w+l ; 
f at=0; 

Sfatl=l-g(jL-l)*conj(g(jL-l)); 
x=((w/2)*(w-5))+4; 
for t=((w/2)*(w-5))+5: ((w/2)*(w-3))+l 
fat=fat+conj (g( ( (N/2) * (N+3) ) -t) ) * . . . 
g(((N/2)*(N+3))-t); 

end 

g( ( (N/2) * (N+3) ) -x) =sqrt (Sf at 1-f at) ; 

else 

continue 

end 

end 

fatp=0.0; 

for iik=N+l:2*N-2 

f atp=f atp+g(iik) *conj (g(iik) ) ; 

end 

a=l-g(l)~2-fatp; 

fatpp=0.0; 

for iiik=2*N:3*N-3 

f atpp=f atpp+g(iiik) *conj (g(iiik) ) ; 

end 

b=l-g(2)~2-fatpp; 

fatig=0.0; 

for pp=2*N:3*N-3 

f atig=f atig+conj (g(pp+(l-N)))*g(pp) ; 

end 

g2N_l_g3N_2=esc_Q (1,2) -f atig; 
c=real(g2N_l_g3N_2) ; 
d=imag(g2N_l_g3N_2) ; 

if d==0.0 

g(2*N-l)=sqrt(a) ; 
g(3*N-2)=sqrt(b) ; 
format short; 

elseif( (c < 1.0e-010) & (c >-0.1e-05)) 
g(3*N-2)=sqrt(b) ; 
g_2N_l_im=-sqrt (a) ; 
g(2*N-l)=complex(0,g_2N_l_im) ; 
format short; 

elseif ((d < 1.0e-010) & (d > -0.1e-05)) 
g(2*N-l)=sqrt(a) ; 
g(3*N-2)=sqrt(b) ; 
format short; 

else 



Uin Polynomial'/ 1 /. 1 /. 1 /. 
AA= [ (4*d~2-4*a*b) * (c~2-a*b) -4*c~2*d~2] ~2 ; 

BB=[(4*d~2-4*a*b)*(c-2-a*b)+4*c~2*d~2]*. . . 
[8*a~2*b*(c~2-a*b-d~2)]-64*c~2*d~2*. . . 
(c~2-a*b)*(2*a~2*b) ; 

CC=[a~2*d~4+a~2*[(c~2-a*b)~2]-2*d~2*a~2*. . . 
(c~2-a*b)]*[(4*d-2-4*a*b)*(c~2-a*b)+. . . 
4*c~2*d~2]+[4*a~2*b*(c~2-a*b-d~2)]~2+. . . 
64*c~2*d~2* (c~2-a*b) * (d~2+a*b) *a~2 ; 

DD=[a~2*d~4+a~2*[(c~2-a*b)-2]-2*a~2*. . . 
d~2*(c~2-a*b)] *4*a~2*b* (c~2-a*b-d~2) ; 

EE=[a~2*d~4+a~2*(c~2-a*b)~2-2*. . . 
a~2*d~2*(c~2-a*b)]~2; 

0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o/o /o /o /o /o/o /o /o /o/o /o /o /o /o/o /o /o /o /o/o /o 

Pii=[AA BB CC DD EE] ; 
rr=roots (Pii) ; 

for ih=l:8 

if isreal(rr(ih)) 
if (rr(ih) < 1.0) & (rr(ih) > 0.0) 
g_2N_l_im=rr (ih) ; 
g3N_2_im=[2*c*g_2N_l_im+. . . 
sqrt((2*c*g_2N_l_im)~2-4*a*. . . 
[b*g_2N_l_im~2+c-2-a*b] )] /(2*a) ; 
g_2N_l_R=-sqrt(a-g_2N_l_im~2) ; 
g3N_2_R=-sqrt(b-g3N_2_im~2) ; 

ccc=g3N_2_R*g_2N_l_R+. . . 

g3N_2_im*g_2N_l_im; 
ddd=g_2N_l_R*g3N_2_im-. . . 

g_2N_l_im*g3N_2_R; 
ai=g_2N_l_R; 
bi=g_2N_l_im; 
g(2*N-l)=complex(ai,bi) ; 
aii=g3N_2_R; 
bii=g3N_2_im; 

g(3*N-2)=complex(aii,bii) ; 
format short; 

else 

continue 

end 

else 
if ih==8 

disp( 'roots=' ) 
disp(rr) 

disp( 'Polynomial roots ') 
disp(' are not real!') 
dispC Check input states!') 
disp ( ' Input states cannot ' ) 
disp(' be linearly dependent.') 
return 
else 

continue 
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end 
end 
end 

end 

dispC ') 

disp( 'Conclusive and inconclusive') 
dispC amplitude probabilities :') 
g 

pause 

'/.'/.'/.'/.Building final vectors Qf V/XVL 
if dim==2*N-l; 

xdim=2*N-l ; 
else 

xdim=dim; 
end 

for uk=l:N; 

for ui=l:xdim; 

Qf (ui,uk)=0.0; 

end 

end 

7, 

jj=o ; 

for iiiv=l:N; 

jj=jj+i; 

for uk=l:N; 

if uk==j j ; 

Qf (uk, iiiv)=g(uk) ; 
else uk<N+l ; 
continue 

end 

end 

end 

S=N; 

uk=0; 

iyy=2*N-l; 
for j=l:N; 
uk=uk+l ; 

if j==l ; 

for zk=S+l:S+(N-j) ; 
Qf (zk,uk)=g(zk) ; 

end 

elseif j<= N-l ; 
xx=uk-2 ; 
for zk=S+l:S+(N-j)+l; 
zmim=S+l ; 

Qf (iyy-(N-uk)+(zk-zmim)-. . . 
xx,uk)=g(zk) ; 

end 

else 

zk=((N/2)*(N+3))-l; 
Qf (N+l,uk)=g(zk) ; 

end 
S=zk; 



end 

Qf 

pause 

y.y.y.y.y.y.y.y.y.y.y.Externai.product a > 1:1:1:1:1:1:1:1:1:11 

A=zeros (xdim) ; 
for ux=l : N 

A=A +Qf ( 1 : xdim , ux) *Q ( 1 : xdim , ux) ' ; 

end 

disp('A=') 
disp(A) 

[V,Sigma,W]=svd(A) ; 

disp('The unitary matrix Ul :') 

U1=V*W 

dispCTest Ul, i.e, (U1*)(U1) = ') 
disp(Ul*Ul') 

disp('The final configuration:') 
disp('[Qlf Q2f Q3f ... QNf ] ' ) 
for iu=l:N 

Qfc(: ,iu)=Ul*Q(l:xdim,iu) ; 
end 

disp(Qf c) 

disp( 'decomposing (Ul)*(Uo)') 

disp('U=Ul*Uo') 

disp(Ul*Uo) 

pause 

y:/X/.Decomposing U=(U1) (Uo Jfflffl! 

U_aux=Ul*Uo; 
k=0; 

for ifg=l:xdim-2 ; 
for j=ifg+l:xdim ; 
k=k+l; 

a=U_aux(ifg,ifg) ; 

b=U_aux(j ,ifg) ; 

c=sqrt (conj (a) *a+conj (b) *b) ; 

a=a/c ; 

b=b/c; 

if ((conj (b)*b)<0. 0001) & (j==ifg+l); 
V=eye(xdim) ; 
R(:,:,ifg,j)=V; 
U_aux=V*U_aux ; 
continue 

elseif ((conj (b)*b)<0. 0001) & (j>ifg+l); 

V=eye(xdim) ; 

V(ifg,ifg)=conj (a) ; 

R(:,:,ifg,j)=V ; 

U_aux=V*U_aux ; 

continue 
else ((conj (b)*b)>0. 00000001) ; 

V=eye(xdim) ; 

V(ifg,ifg)=conj (a) ; 

V(ifg, j)=conj (b) ; 

V(j,ifg)=b; 

V(j,j)=-a; 
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R(: , : ,ifg,j)=V ; 
U_aux=V*U_aux ; 
continue 

end 
end 
end 

V_aux=U_aux; 
if g=xdim-l ; 

j=ifg+i; 

k=k+l; 
V=eye(xdim) ; 

V(ifg,ifg)=conj (V_aux(ifg,ifg)) ; 
V(ifg, j)=conj (V_aux(j ,ifg)) ; 
V(j ,ifg)=conj (V_aux(ifg, j)) ; 
V(j , j)=conj (V_aux(j , j)) ; 
R(:,:,ifg,j)=V; 
for ifg=l:xdim-l ; 

for j=if g+1 : xdim ; 



if j==ifg ; 
continue 

else 

Rotation(: , : ,ifg,j)=R(: , : ,ifg, j) ' ; 
disp( 'Rotation' ) ,disp( [if g j]) 
disp(Rotation( : , : , if g, j ) ) 
end 

end 
pause 

end 

disp('Test Rotations') 

disp('R_(d-l,d)R_(d-2,d) . . .R(1,3)R_(1,2)U=') 

U_f inal=V*V_aux; 

disp(U_f inal) 
llll END '/.'/.'/.'/.'/.'/.'/. 



[1] A. Peres, Quantum Theory: Concepts and Methods 
(Kluwer Academic Publishers, Dordrecht, 1993). 

[2] C. H. Bennett, Phys. Rev Lett. A 68, 3121 (1992). 

[3] I. D. Ivanovic, Phys. Lett. A 123, 257 (1987). 

[4] D. Dieks, Phys. Lett. A 126, 303 (1988). 

[5] A. Peres, Phys. Lett. A 128, 19 (1988). 

[6] G. Jaeger, A. Shimony, Phys. Lett. A 197, 83 (1995). 

[7] A. Chefles, Phys. Lett. A 239, 339 (1998). 

[8] Y. C. Eldar, IEEE Trans. Inform. Theory, 49, 446- 
456 (2003). 

[9] Sonja Franke- Arnold, E. Andersson, S. M. Barnett, 
S. Stenholm, Phys. Rev. A 63, 052301 (2001). 

[10] Y. Sun, M. Hillery, J. A. Bergou, Phys. Rev. A 64 
022311 (2001). 

[11] L. Roa, J. C. Retamal, C. Saavedra, Phys. Rev. A 
66, 012103 (2002). 



12] S. Boyd, L. Vandenbergue, "Semidefinite Program- 
ming", Siam Review, March 1996. 

13] M. A. Neumark, Izv. Akad. Nauk. SSSR, Ser. Mat. 
4, 53 277 (1940). 

14] R. A. Horn, C. Johnson, Matrix Analysis, (Cambridge 
University Press, 1991). 

15] M. Reck, A. Zeilinger, et. al., Phys. Rev. Lett. 73 , 
58 (1994). 

16] M. Dusek, M. Jahma, N. Liitkenhaus, Phys. Rev. A 
62, 022306 (2000). 

17] S. J van Enk, Phys. Rev. A 66, 042313 (2002). 

18] C.H. Bennett and G. Brassard, in Proceedings of the 
IEEE International Conference on Computer, Sys- 
tems, and Signal Processing, Bangalore, India (IEEE, 
New York, 1984), pp. 175-179. 



